Stability of Repulsive Bose-Einstein Condensates in a Periodic Potential 
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The cubic nonlinear Schrodinger equation with repulsive nonlinearity and an elliptic function 
potential models a quasi-one-dimensional repulsive dilute gas Bose-Einstein condensate trapped in 
a standing light wave. New families of stationary solutions are presented. Some of these solutions 
have neither an analog in the linear Schrodinger equation nor in the integrable nonlinear Schrodinger 
equation. Their stability is examined using analytic and numerical methods. All trivial-phase stable 
solutions are deformations of the ground state of the linear Schrodinger equation. Our results show 
that a large number of condensed atoms is sufficient to form a stable, periodic condensate. Physically, 
this implies stability of states near the Thomas-Fermi limit. 



I. INTRODUCTION 



Recent experiments on dilute-gas Bose-Einstein con- 
densates (BECs) have generated great interest in macro- 
scopic quantum phenomena in both the theoretical 
and experimental physics community. Such BECs are ex- 
perimentally realized when certain gases are super-cooled 
below a critical temperature and trapped in electromag- 
netic fields j^. Many BEC experiments use harmonic 
confinement. Recently, however, there has been much 
interest in sinusoidal confinement of repulsive BECs us- 
ing standing light waves. Such BECs have been used 
to study phase coherence |^-|^ and matter-wave diffrac- 
tion . They have also been predicted to apply to quan- 
tum logic , matter- wave transport [|lO| , and matter- 
wave gratings. In this paper, we consider the dynamics 
and stability of repulsive BECs trapped in standing light 
waves. 

A mean-field description for the macroscopic BEC 
wavefunction is constructed using the Hartree-Fock ap- 
proximation [ prf and results in the Gross-Pitaevskii equa- 
tion The dimensions of the BEC play an impor- 
tant role: ID, 2D, and 3D BECs all behave in a radi- 
cally different manner In the quasi- ID regime, 
the Gross-Pitaevskii equation reduces to the ID nonlin- 
ear Schrodinger equation (NLS) with an external poten- 
tial. This regime holds when the transverse dimensions 
of the condensate are on the order of its healing length 
and the longitudinal dimension is much longer than its 
transverse dimensions p6| , p7[ |. In this regime the BEC 
remains phase coherent and the governing equations are 
one-dimensional. This is in contrast to a truly ID mean- 
field theory which requires transverse dimensions on the 
order of or less than the atomic interaction length. 

The recent trapping of a BEC in a hollow blue-detuned 
laser beam demonstrates that a quasi-lD BEC is 
experimentally realizable. A variety of other experi- 



ments |]jT8[|T|,|2|-||l are also modeled by the ID NLS 
with an external potential. Upon rescaling, the govern- 
ing evolution is given by 



(1) 



where ip{x,t) is the macroscopic wave function of the con- 
densate and V{x) is an experimentally generated macro- 
scopic potential. Confinement in a standing light wave re- 
sults in V{x) being periodic. In a recent experiment 
a shallow harmonic potential was applied in addition to 
a standing light wave. The standing light wave in this 
case was sufficiently intense so that the condensate was 
strongly localized in each well. This is referred to as the 
tight-binding regime. Additionally, the apparatus was 
tilted vertically so that gravity caused tunneling between 
wells. Our theoretical findings consider the complimen- 
tary experiment in which the condensate is free to move 
between wells. With the advent of quasi-lD, cylindrical 
geometries | |2^ ], additional harmonic confinement is no 
longer necessary and the BEC dynamics considered here 
are applicable. 

To model the quasi-lD confinement produced by a 
standing light wave, we use the periodic potential 



V{x) = -Va sn^{x,k) 



(2) 



where sn(x, k) denotes the Jacobian elliptic sine func- 
tion 1^ with elliptic modulus < fc < 1. In the limit 
k = the potential is sinusoidal and thus V{x) is a stand- 
ing light wave. For intermediate values (e.g. k < 0.9) 
the potential closely resembles the sinusoidal behavior 
and thus provides a good approximation to a standing 
light wave. Finally, for k — 1~, V{x) becomes an ar- 
ray of well-separated hyperbolic secant potential barriers 
or wells. The potential is plotted in Fig. |l| for values of 
fc = 0, 0.9, 0.999 and 0.999999. Only for fc very near one 
(e.g. fc > 0.999) does the solution appear visibly elliptic. 
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k=1 




4>{x,t) = r{x) eKp(—iLut + i9(x)). 



(3) 



2K(k)x/7r 

FIG. 1. The sn^{x,k) structure of the potential for 
varying values of k. Note that the a;-coordinate has been 
scaled by the period of the elliptic function. This period 
approaches infinity as fc — > 1. Since sn{x,k) is periodic in 
X with period 4:K{k) = 4 J^^^ da/\/ 1 — fc^ sin^ a, V{x) is 
periodic in x with period 2K{k). 



The freedom in choosing k allows us to consider 
much more general potentials than considered previ- 
ously [psf-^ and allows for great flexibility in consid- 
ering a wide variety of physically realizable potentials. 

The paper is outlined as follows: in the next section we 
derive and consider various properties and limits of two 
types of explicit solutions of Eq. (|^) with (||) . Section III 
develops the analytic framework for the linear stability 
properties of the new solutions of Sec. II. The stability 
results are confirmed by numerical computations. In cer- 
tain cases, the stability analysis is intractable and we rely 
solely on simulations to determine stability. We conclude 
the paper in Sec. IV with a brief summary and highlights 
of the primary results of the paper and their consequences 
for BEC dynamics and confinement. 



II. STATIONARY SOLUTIONS 

Equation (0) with V{ x) = is an integrable equa- 
tion and many explicit solutions corresponding to vari- 
ous boundary conditions arc known. A comprehensive 
overview of these solutions is found in If V{x) ^ 0, 
the NLS is not integrable. In this case, only small classes 
of explicit solutions can most likely be obtained. Our 
choice of potential (j^) is motivated by the form of the 
stationary solution of the NLS with V{x) = 0. An 
overview of these stationary solutions and their proper- 
ties is found in p6[ |. At present, we restrict our attention 
to stationary solutions of Eq. (Q), i.e., solutions whose 
time-dependence is restricted to 



If Ox = 0, then the solution is referred to as having trivial 
phase and we choose 0{x) = 0. Substituting the ansatz 
Eq. (H) in Eq. (|^) and dividing out the exponential factor 
results in two equations: one from the real part and one 
from the imaginary part. The second equation can be 
integrated: 



{x) = c 



dx' 

ipl ( nnl 



(4) 



where c is a constant of integration. Note that d(x) is a 
monotonous function of x. Substitution of this result in 
the remaining equation gives 



{x) = — 



c r{x)r"[x) 



(5) 



The following subsections describe two classes of solu- 
tions of this equation. 



Type A 

1. Derivation 

For these solutions, r'^{x) is a quadratic function of 
sn(a;, k): 



r'^{x) = A sn'^{x,k) + B. 



(6) 



Substituting this ansatz in Eq. (|5|) and equating the co- 
efficients of equal powers of sn(a;, k) results in relations 
among the solution parameters uj,c,A and B and the 
equation parameters Vq and k. These are 



(7a) 



c' = B 1 



B 



Vo + k^ 



A = Vo + k'^. 



{Vo + P + Bp) , (7b) 
(7c) 



For a given potential V{x), this solution class has one 
free parameter B which plays the role of a constant back- 
ground level or offset. The freedom in choosing the po- 
tential gives a total of three free parameters: Vq, k and 
B. 

The requirements that both r'^{x) and are positive 
imposes conditions on the domain of these parameters: 



Vq > -P, B>0, or 



(8a) 



Vq < -k\ -{VQ + k')<B<-{l + -^ \ . (8b) 



Vq 



2 



B=-(Vo+k2)/k^ . B 





Vn 



-k^ 



FIG. 2. The region of validity of the solutions of Type 
A is displayed shaded for a fixed value of k. The edges of 
these regions correspond to various trivial phase solutions. 



The region of validity of these solutions is displayed in 
Fig. |. 

For typical values of Vq, A: and B, the above equations 
give rise to solutions of Eq. which are not periodic in 
x: r{x) is periodic with period 2K{k), whereas exp(i0(x)) 
is periodic with period T = f?~^(27r). In general these 
two periods 2K{k) and T are not commensurable. Thus, 
requiring periodic solutions results in another condition, 
namely 2K{k) /T — p/q, for two positive integers p and q. 
The most convenient way to express this phase quantiza- 
tion condition is to assume the potential (i.e., Vq and k) is 
given, and to consider values of B for which the quantiza- 
tion condition is satisfied. Introducing (3 — B/ (Vb + fc^), 
we find 



± 



dx 



sn(a;, fc)^ + (3 q 



(9) 



This equation is solved for /3, after which B = /3(Vo-f fc^). 
For numerical simulations, the number of periods of the 
potential is set. This determines q, limiting the number 
of solutions of Eq. . Solutions with the same period- 
icity as the potential require p/q — 1. 

Note that solutions of Tvpe A reduce to stationary so- 
lutions of Eqs. and (||) with Vq = 0. Furthermore, 
all stationary solutions of the integrable equation are ob- 
tained as limits of solutions of Type A. 



2. Limits and Properties 

The properties of these solutions are best understood 
by considering their various limit cases. 

The trivial phase case: The solutions of Type A 
have trivial phase when c = 0. Since has three factors 
which are linear in B (see Eq. (^), there are three choices 
of B for which this occurs: B = 0, B = — (Vq -I- k'^) and 
B — — (Vq -|- fc^)/fc^. These possibilities are three of the 
four boundary lines of the region of validity in Fig. |^. 
Note that the remaining boundary line (Vb = — fc^) cor- 
responds to r^(x) = B, which gives rise to a plane wave 



solution. Using Jacobian elliptic function identities [g7j , 
one finds that the three other boundary lines give rise to 
simplified solution forms: B = Q gives 



(10) 



r{x) = \/Vq + k^ sn(x, k), uj ^ 

B = -{Vo + P) gives 

r(x) = y^O^T+fc^y cn(a;,fc), uj^'^-Vo-k', (11) 

where cn(a::, k) denotes the Jacobian elliptic cosine hmc- 
tion. Lastly, B — —{Vo + k^)/k'^ gives 



(12) 

where dn(a;, k) denotes the third Jacobian elliptic func- 
tion. Solution ( p^ ) is valid for Vb > — fc^, whereas the 
other two solutions (11) and ( p^ ) are valid for Vb < —k^. 
The amplitude of these solutions as a function of poten- 
tial strength Vb is shown in Fig. 0. 
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FIG. 3. The amplitude of the trivial phase solutions of 
Type A versus the potential strength Vb. 
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f(x)=-(l+V/cT dn (x,k) 



r (x)=-(V +kl cn"(x.k) 



V(x) - -V, 




2K{k) 
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FIG. 4. Trivial phase solutions for k = 0.5. V{x) is 
indicated with a solid line. For the top figure Vb = 1. For 
the bottom figure Vo = —1. 
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(a) V„=1,B=1 



(b) V„=1,B=0.1 




It 

FIG. 5. Phase and amplitude of the trigonometric solu- 
tions. For all these figures, the solid line denotes V{x), the 
dashed line is r{x) and the dotted line is 9{x)/{2tt). Note 
that 9{x) becomes piecewise constant, as B approaches the 
boundary of the region of validity. Far away from this 
boundary, 0{x) is essentially linear. 

Both cn(a;, k) and sn(a;, k) have zero average as func- 
tions of X and he in [—1,1]. On the other hand, dn(a;, k) 
has nonzero average. Its range is — k'^, 1] . Further- 
more, cn(x, k) and sn(x, k) are periodic in x with period 
4K{k), whereas dn{x,k) is periodic with period 2K{k). 
These properties matter greatly for the stability analysis, 
as will be seen in Section 3. Some solutions with trivial 
phase are shown in Fig. ^ 

The trigonometric limit: In the limit fc ^ 0, the 
elliptic functions reduce to trigonometric functions and 
V{x) = -Va sin^{x) = (Vb/2) cos(2a:) - Vo/2. Then 

r^{x)^Vo s\u^{x) + B, uj=\ + B. (13) 



In this case, the phase integral Eq. (H) results in 



tan(6l(x)) ±yjl + VqIB tan(a;). 



(14) 



Note that this formula guarantees that the resulting so- 
lution is periodic with the same period as the potential, 
so no phase quantization is required. In the trigonomet- 
ric limit, the wedge between the two regions of validity 
in Fig. ^disappears. This is no surprise, as in this limit, 
dn(a:, /c) — *■ 1, and the third trivial phase solution reduces 
to a plane wave solution. The cornerpoint of the region 
of validity also moves to the origin. Some trigonometric 
solutions are illustrated in Fig. ^j. 

The solitary v^rave limit: k — \. In this limit the 
elliptic functions reduce to hyperbolic functions. Specifi- 
cally, sn(x, /c) = tanh(2;). Hence in this limit, the poten- 
tial has only a single well or a single peak. Then Vb < 



gives rise to a repulsive potential, whereas Vq > gives 
rise to an attractive potential: V{x) = — Vq tanh^(a;). 
In this case the phase Q{x) of Eq. (Q) can be calculated 
explicitly: 



r^(x) = (Vb + 1) tanh2(a;) + B, 



(15a) 

0{x) — ^/ -j-;-^^— j-a;-|-arctan^Y^^^^tanh(a;)^ , (15b) 



B 



Vb + 1 



which is valid for Vb > — 1 and B > 0. This solution is a 
stationary solitary wave of depression on a positive back- 
ground. It is a deformation of the gray soliton solution 
of the NLS equation with repulsive nonlinearity. Note 
that these solutions can exist with an attractive poten- 
tial provided — 1 < Vq < 0. Two solutions with repulsive 
potential are illustrated in Fig. |6^-b. Another solution is 
valid for S = -(Vb -I- 1) > 0: r{x) = ^-(Vb -h 1) sech(x) 
and Q(x) — 0. This solution represents a stationary ele- 
vated solitary wave. It is reminiscent of the bright soli- 
ton solution of the NLS equation with attractive nonlin- 
earity. This solution is shown in Fig. A surprising 
consequence of considering Eq. (|]) is that the potential 
strength Vb acts as a switch between the equation with re- 
pulsive and attractive nonlinearity, as illustrated by these 
solitary wave solutions. 

Understanding the solitary wave limit facilitates the 
understanding of what occurs for fc 1. In this case 
the solutions of Type A reduce to a periodic train of soli- 
tons with exponentially small interactions as illustrated 
in Fig. |d. 



(a) V„=1, B=1 



(b) V„=1,B=0.1 





5 -5 



c) V„=-2, B=1 (d) V„=1 , 8=2, k=1 -1 0" 



1.5 - 




5 -4K(k) 



4K(k) 



FIG. 6. Solutions with k=l (a,b,c) or fc ^ 1 (d). The 
solid line denotes Vix), the dashed line is r(x) and the dot- 
ted line is 9{x)/2-k. In (d), a value of fc = 1 — 10^^^ was 
used. 
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Type B 

1. Derivation 

For these solutions, r'^{x) is linear in cn(a:;, fc) or 
dn{x,k). First we discuss the solution with cn(x, A;). The 
quantities associated with this solution will be denoted 
with a subindex 1. The quantities associated with the 
dn(a;, k) solution receive a subindex 2. 

Substituting 



rl{x) — ai cn{x, k) + bi 



(16) 



in Eq. (^) and equating different powers of cn(a;, k) gives 
the relations: 



(17a) 
(17b) 



zj = ^il6al~k^){16al + e^k^), (17c) 



bi 



4a? 
/fc2 



(17d) 



The class of potentials Eq. (||) is restricted by the first 
of these relations so that Vq is in the narrow range 
< Vb < 0. The solution class now depends on 
one free amplitude parameter ai and the free equation 
parameter fc. 

The region of validity of this solution is, as before, de- 
termined by the requirements c? > and rKx) > 0: 



fc2 



(18) 



The class of potentials (||) is restricted as for the previous 
solution by the first of these relations. The solution class 
again depends on one free amplitude parameter 02 and 
the free equation parameter k. 

The region of validity of this solution is once more de- 
termined by the requirements c| > and r|(x) > 0: 



, , 1 Vl - 
laal > - or < 02 < . 



(22) 



The period of r2{x) is equal to the period of the po- 
tential. Requiring periodicity in x of this second solution 
of Type B gives 



± 



^/{16al-l){16al + k^~l) 



K(k) 



dx 



P 



4a2-t-dn(a;, fc) 



(23) 



For given fc and integers p, q, this equation needs to be 
solved to determine 02- 

In contrast to solutions of Type A, solutions of Type 
B do not have a nontrivial trigonometric limit. In fact, 
for solutions of Type B, this limit is identical to the limit 
in which the potential strength Vb — — 3fc^/8 approaches 
zero. Thus it is clear that the solutions of Type B have 
no analogue in the integrable nonlinear Schrodinger equa- 
tion. However, other interesting limits do exist. 



2. Limits and Properties 



The period of ri (x) is twice the period of the potential. 
Requiring periodicity in x of this first solution of Type B 
gives 



± 



vW^^Fycsf+T^P) /■'^^''^ dx 



An 



P 



4:(3i + k cn{x, fc) q 



(19) 



For given fc and integers p, q, this equation is solved for 
from which oi = f3ik/A. 
The dn(a::, fc) solutions are found by substituting 



r\{x) = 02 dn(a::, fc) + 62 



(20) 



in Eq. (H). Equating different powers of dn(a;, fc) imposes 
the following constraints on the parameters: 

3,, 



0J2 = -(l + fc2) + 6a2, 



(21a) 
(21b) 



■-{\&al-l){\&al + k'^ -I) , (21c) 



62 = 4a2. 



(21d) 



The trivial phase case: Trivial phase corresponds 
to c = 0. This occurs precisely at the boundaries of the 
regions of validity. For the first solution of Type B, there 
are two possibilities: ai = k"^ /A or ai = —k"^ /A. By 
replacing a; by a; -I- 2K{k), one sees that these two pos- 
sibilities are completely equivalent, so only the first one 
needs to be considered: 



fc2 

^iW = ^(l + cn(x,fc)) 



1 



fc2 

2 ■ 



(24) 



For the second solution, there are four possibilities: 02 = 
1/4, a2 = -1/4, 02 = and 02 = Vl - k'^/A. The third 
one of these results in a zero solution. The others give 
interesting trivial phase solutions. For 02 = 1/4, 



The case 02 = —1/4 gives 



1 k' 



r2(a;) = -(1 - dn(a;, fc)) , a;2 = - -I- 



(26) 
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Finally, for 02 = \/r~~P/4, 



(dn(x,/c) + \/l-fc2), 



W2 = 



(27) 



0.5 



(2) : 






(4) 


















(3) 



2K(k) 



4K(k) 




4K(k) 



2K(k) 



4K(k) 



FIG. 7. Solutions of Type B with trivial phase. The 
figures correspond to, from top to bottom, k — 0.5, k = 0.9 
and k = 0.999. The potential is indicated with a solid line. 
The other curves are: (1) |ri(a;)| with 02 = (2) r2{x) 

with 02 = 1 /4, (3) |r2(a;)| with 02 = -1/4 and (4) r2{x) 
with a2 = Vl — fc^/4. 
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5 -2.5 


2.5 E 



FIG. 8. Solitary wave solutions of Type B. The poten- 
tial is indicated with a solid line. The dashed line solution 
corresponds to a = 0.3, the dotted line to a = —0.3. 

These solutions are shown in Fig. |^. 

The solitary wave limit: In this limit Vq — —3/8 
and the potential is V{x) — 3tanh(x)^/8. Since both 
cn(a:, fc) = sech(a:) and dn(x, fc) — sech(x) when k = 1, 
ri(x) = r2{x) in the solitary wave limit. Their ranges 
of validity also share the same limit: \a\ > 1/4 with 
a = ai =02. The phase can be calculated explicitly and 
the solitary wave solution of Type B is 



0ix) 



1 / /4a-l x\ , , , 

— =Farctan I W ^ tanh — 1 . (28bj 



The region of validity consists of two separated regions: 
a > 1/4 and a < —1/4. In the first region, the solution 
is a stationary elevated solitary wave with a constant 
background da^ . In the second region a < — 1/4 and the 
solution is a stationary solitary wave of depression with a 
constant background 4a^. These solitary wave solutions 
are illustrated in Fig. |^ 

As for the solutions of type A, the solitary wave limit 
gives an idea of the behavior of the solution for values of 
k —^ 1, where the solution behaves as a periodic array of 
solitons with exponentially small interactions. 



III. STABILITY 

We have found a large number of new solutions to the 
governing Eqs. (|l]) and (||). However, only solutions that 
are stable can be observed in experiments. In this sec- 
tion, we consider the stability of the different solutions. 
Both analytical and numerical results are presented for 
the solutions with trivial phase. In contrast, only numer- 
ical results are presented for the nontrivial phase cases. 

We first consider the linear stability of the solution (||). 
To do so, consider perturbations of the exact solutions of 
the form 

ip{x, t) = [r{x) + e(j)[x, t)) e-x:p[i{e{x) - uot)] (29) 

where e <C 1 is a small parameter. Collecting terms at 
0{e) gives the linearized equation. In terms of the real 
and imaginary parts U — {Ui, U2Y — (i?e[0], /to[0])* the 
linearized evolution is given by: 



JL\J = J 



-S 



U, 



(30) 



where 
L+ = 
L_ = 
S = 

and J = 



r'^{x) 



2 V r^x) 

— 5x — 

r{x) r{x) ' 

-1 

1 



+ 3r^{x)+V{x)-uj, (31a) 
+ r^{x) + V{x) - LO, (31b) 
(31c) 



is a skew-symmetric matrix. The 

operators L, and are Hermitian while S is anti- 
Hermitian. Considering solutions of the form \J{x,t) = 
\J{x) exp(At) gives the eigenvalue problem 



r^(x) = 4a^ -I- a sech(a;) 



(28a) 



£U = AU, 



(32) 
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where £ = JL and A is complex. If all A are imaginary, 
then linear stability is established. In contrast, if there 
is at least one eigenvalue with a positive real part, then 



instability results. Using the phase invariance V' 
of Eq. (0), Noether's theorem Q gives 





r{x) 



= 0, 



(33) 



which implies that L-r{x) — 0. Thus A = is in the 
spectrum of L_. For general solutions of the form (|^), 
determining the spectrum of the associated linearized 
eigenvalue problem ( ^o|) is beyond the scope of current 
methods. However, some cases of trivial phase solutions 
(c = 0) are amenable to analysis. 

The Hermitian operators L± are periodic Schrodinger 
operators and thus the spectra of these operators is real 
and consists of bands of continuous spectrum contained 
in [A±,cx3) Here A± denote the ground state eigen- 
values of L± respectively. They are given by 



A+ = inf 



(34) 



where \\4>\\'^ — From the relation — L_ + 2r^(x) 

it follows that A+ > A_. Also A_ < since A = is an 
eigenvalue of L_. 

If A+ > 0, then is positive, so we can define the 

positive square root, L^, via the spectral theorem [Q, 

and hence the Hermitian operator H = L^L^L^ can be 
constructed. The eigenvalue problem for C in Eq. (^2|) is 
then equivalent to 



(iJ + A2)(^ = 0, 



(35) 



with ip = L'^Ui. Denote the left-most point of the spec- 
trum of H by fiQ. If /iQ > then A^ < and the eigenval- 
ues of C are imaginary and linear stability results. Since 

11 1 

H = L^L^L^ and is positive, ^j-q > ii and only if 
L_ is non-negative. In contrast, if /Zq < then A^ > 
and £ has at least one pair of real eigenvalues with op- 
posite sign. This shows the existence of a growing mode 
leading to instability of the solution. 

Three distinct cases are possible for linear stability 

• If r{x) > then r{x) is the ground state of L_ [^ , 
and Eq. (|3^ ) implies A_ = and hence A+ > 0. 
Thus the solution (||) is linearly stable. 

• If r{x) has a zero, it is no longer the ground 
state Q and A_ < 0. Thus there exists a tpo such 

that (tpolL-lipo) < 0. If in addition A+ > 0, then 

_i _i 

we can construct (f> — Lj^^ipo/\\I^+'^4'o\\ which gives 
{(f>\H\(j)) < 0. Hence /io < and £ has positive, real 
eigenvalues so that the solution (|^) is linearly un- 
stable. 



• For A_ and A+ both negative the situation is indef- 
inite and our methods are insufhcient to determine 
linear stability or instability. 

In what follows, these results are applied to the Type A 
and B trivial phase solutions constructed in the preceding 
section. Specifically, we construct the operators L_ 

T Im(X) 



dn(x,k) 



cn(x,k) 
— i 



sn(x,k) 



-1/2 



(k-l)/2 



FIG. 9. The spectrum of 
trivial phase solution. 



for the Type A cn(a;, k) 



and for each solution, which allows us to use one 
of the above criteria. The analytical results are accom- 
panied by direct computations on the nonlinear govern- 
ing Eqs. (|l|) and (||). For all computational simulations, 
twelve spatial periods are used. However, to better il- 
lustrate the dynamics, typically four spatial periods are 
plotted. Moreover, all computations are performed with 
white noise included in the initial data. 



A. Trivial Phase: Type A 



1. cn(x,k) 



For the cn(x, k) solution the L± operators are 



- {2Vo + 3fc2)cn2(x, A:) + fc^ - 1 (36a) 



2 



fc^cn^(a;, k) + k^ 



(36b) 



with Vq < —k^. Note that L_, which is independent of 
Vb, is the classical 1-gap Lame operator The spec- 

trum of L_ can be calculated explicitly. The ground state 
eigenvalue is A_ = {k"^ — l)/2 with associated eigenfunc- 
tion dn(x, k). The elliptic functions cn(a:;, k) and sn(a;, k) 
are also eigenfunctions of They are the first and sec- 
ond excited state and have eigenvalue and k"^ /2 respec- 
tively. These are the only eigenvalues and the spectrum 
consists of the bands [{k'^ ~ l)/2,0] U [fc2/2,oo). The 
spectrum is illustrated in Fig. |^. 

Since dn(x, k) > Q and A_ = (/c^ - l)/2 < the argu- 
ments of the previous section imply that the cn(x, k) wave 
is unstable whenever the operator i+ > 0. It is clear 
from Eq. (|6|a) that i+ is positive if Vq < -{k"^ + 1/4) 
and fc^ > 1/2. Thus, the cn(a;, fc) wave is unstable for 
parameter values in this region. Moreover, this region 
can be enlarged to Vb < -{k^ + l)/2 and k^ > 1/2 by 
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observing that the ground state eigenvalue of an opera- 
tor Lq + 7L1 is a convex function of 7: A = A(7). This 
follows from the fact that the ground state eigenvalue is 
the minimizer of the Rayleigh quotient. Let a G [0, 1], 
then 




FIG. 10. Unstable evolution of a Type A cn(a;, k) solu- 
tion given by Eq. ( pA[ ) over 40 time units with k = 0.5 and 
Vo = -1.0. 




I I 

-1 1 wavenumber 

FIG. 11. Wavenumber spectrum evolution of a Type A 
cn{x, k) solution given by Eq. (|l^) over 100 time units with 
k = 0.5 and Vo = —0.55. The modal evolution shows the 
band of unstable modes which result from starting with the 
unstable cn(a;, k) solution. This shows that the instability 
occurs in a neighborhood of the dominant wavenumber of 
the stationary solution. 

A (a7i -I- (1 - 0)^2) 

= inf ((/)|a(Lo + 7iii) + (1 - a)(^o +72^1)10) 
11*11=1 

> a inf {(f)\LQ+^iLi\(f)) + (1 — a) inf {(j)\LQ+^2Li\(j)) 
11011=1 Wll=i 

= aA(7i)-f (l-a)A(72). (37) 



Now consider the ground state eigenvalue A_|_ = A_|_(Vb) 
and note that A+(-fc2) = {p _ l)/2 and A+(-3A;2/2) = 
fc^ — 1/2. The line through these two points is given by 
A+(Vo) = -Vo - (1 + A:2)/2, so by convexity A+(Vb) > 
-Vo -{1 + P)/2 for Vo e h3fcV2,-fc2]. Thus, A+ = 
MVo) > if < -il + P)/2. 

If k'^ < 1/2, less is known. However the results of 
Weinstein and Keller [|4j show that the ground state 
eigenvalue grows as A+ « (2(1 - A:^)|V"o|)^/^ + fc^ - 1/2 
for — Vo ^ 1. Hence for fc^ < 1/2 instability occurs for 
sufficiently negative Vq. 

The most unstable modes of the cn(a;, k) solution can 
be determined perturbatively when e = — 2(Vb-l-fc^) <C 1. 
This corresponds to a solution with small amplitude. 
Since = + ecn'^{x,k), it follows that L+ is not 
necessarily positive, disallowing the construction of H in 
Eq. (H). However, from Eq. (H), L+L_U2 = -A^C/z, 
which offers an alternative to Eq. ( ^5|) to calculate the 
spectrum of Eq. (|3^). Let A = iu+eXi and U2 = (j>u+£(t)i, 
where ly is an eigenvalue of L_ and 4>i, is its associated 
normalized eigenfunction. Then a first order calculation 
using time-independent perturbation theory gives 

= -ly^ ~ siy{4>^\cn'^{x,k)\(l}^) . (38) 

Thus, A^ > only if -e {(j)iy\cn^{x, k)\(j),,) < ly < 0. 
Hence, only modes </)^ with v in this band near zero are 
unstable. For these unstable modes, the eigenfunction 
(pi, is approximately the zero mode cn(a;,fc). Thus the 
onset of instability in the Fourier domain occurs near the 
wavenumbers of the cn(x,k) solution. This is character- 
istic of a modulational instability 

To illustrate this instability, we display in Fig. |lO| the 
evolution of a cn(2:, k) solution over the time interval 
t e [0,40] for Vo = -1.0 and k = 0.5. The solution goes 
quickly unstable with the instability generated near the 
first wavenumber. This agrees with the analytical predic- 
tion. It is illustrated in the evolution of the wavenumber 
spectrum in Fig. |ll|. Here a close-up of the spectrum 
near wavenumber one is shown. This shows that the in- 
stability indeed occurs in a neighborhood of the dominant 
wavenumber of the stationary solution. 

2. sn{x, k) 

For the sn(x, k) solutions the L± operators are given 

by 

L+ = -^5,, + (3fc2 + 2VoW{x, k) - (39a) 
i- = -^a,, + fc2sn2(x,fc)-i^. (39b) 

Again is a 1-gap Lame operator, differing from for 
the cn(x, k) solution only by a constant. The spectrum 
is given by [— /c^/2, — 1/2] U [0, 00). It again follows from 
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the work of Weinstein and Keller that for sufRciently 
large values of Vq the ground state eigenvalue of L_ is 
approximately given by 



B. Trivial Phase: Type B 



A+ « (2^0)^ - 



(40) 




1. cn{x,k) 



FIG. 12. Stable evolution of the Type A dn(a;, k) solu- 
tions given by Eq. ( p^ ) over 80 time units with k = 0.5 and 
Vo = -1.0. 



and thus L_|_ is positive definite for sufficiently large Vq. 
This, in turn, implies instability of the sii[x, k) solution 
for sufficiently large Vb, which corresponds to large am- 
plitude solutions. The sn(a;, k) solution goes quickly un- 
stable in a similar fashion to the cn(x, fc) solution (see 

Fig. m. 



The Type B trivial phase solution is obtained for 
Oi = ±fc^/4 and corresponding amplitude \r{x)\ = 
(k/2)^Jl + cn(x, k). The solution r{x) is not strictly pos- 
itive. The operator L+ is 



1 3 

- + -fc^sn^(x, A:)-)-3aicn(a;,fc). (41) 
8 8 



3. dn{x,k) 



From the previously established results, linear stabil- 
ity for the dn(a;, k) solutions follows immediately since 
r(x) > 0, because dn(a;, k) has no zeros. Thus in con- 
trast to the cn(a;, k) and sn(x, k) solutions, the dn(a;, fc) 
solutions given by Eq. (12) are linearly stable. Figure 
displays the evolution of a dn(a;, fc) solution over the time 
interval t G [0,80] for Vq = -1.0 and fc = 0.5. Although 
noise was added to the initial data, the solution shape 
persists and the solution is stable, as predicted analyt- 
ically. For this case, the wavenumber spectrum is sup- 
ported primarily by three modes: the zero mode which 
determines the offset, and two other modes which de- 
termine the oscillation frequency of the dn(x, fc) solu- 
tion. Even with large perturbations, this solution per- 
sists. This indicates that the offset of a solution is im- 
portant for its stability. This observation is reconfirmed 
for other stable solutions below. 



Thus we find that the situation is indeterminate. 



Numerical simulations for the Type B cn(a:, fc) solu- 
tions given by Eq. ( |2^ ) are illustrated in Fig. |l^. This 
figure displays the evolution of the cn(a;, fc) branch of so- 
lution for fc = 0.5 (top panel) and fc = 0.999 (bottom 
panel) over the time interval t e [0, 800] and t e [0, 400] 
respectively for Vq — — 3fc^/8. For both fc = 0.5 and 
fc = 0.999 the solutions are unstable, but this instabil- 
ity manifests itself only after several hundred time units. 
Figure ^ shows the evolution of the wavenumber spec- 
trum for both these cases. For fc = 0.5, the onset of 
instability occurs near wavenumber one as is the case of 
Type A solutions. After 800 time units, the wavenumbers 
have only just begun to spread, causing the solution to 
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X 27 




-400 



-36 



FIG. 13. Unstable evolution of the Type B cn(x, k) so- 
lutions given by Eq. for k = 0.5 (top panel) and 
k = 0.999 (bottom panel) for ai = k^/A and Vo = -3k^/8. 



800 
t 



400 r 
t 



mode 



mode 



FIG. 14. Wavenumber spectrum evolution of a Type B 
cn{x,k) solution given by Eq. ^) for Oi = P/4 and cor- 
responding to A: = 0.5 (left panel) and k = 0.999 (right 
panel) of Fig |l3| The evolution shows that the unstable 
band of modes is generated near wavenumber one and for 
k = 0.999 near wavenumber one and its harmonics. 




FIG. 15. Stable evolution of a Type B dn(a;, k) solution 
given by Eq. for k = 0.999 and 02 = 1/4. 



destabilize. For k — 0.999, the solution is composed of 
a much larger number of wavenumbers which destabi- 
lize much more quickly than the k = 0.5 case. Here 
the instability is generated near wavenumber one and its 
harmonics. 



2. dn{x,k) 

The trivial phase dn(a;, k) solution requires c = which 
is achieved for 02 = ±1/4, 02 = 0, or 02 = Vl — k'^/A. 
Thus three distinct parameter regimes need to be consid- 
ered. The relevant operators in this case are 



1 fc^+1 3fc^ 

-dl — + 6al+—sn\x,k) + 3a2dn{x,k), (42a) 

2 o o 

1 + 1 '\k'^ 

- -dl 2a\v— sn^{x,k)+a2 dn {x,k). (42b) 

2 8 8 



The case 02 = 1/4 gives L^r{x) — with r{x) > 0. 
Hence from the linear stability criteria, these waves are 
stable for all values of k. As with the a2 = 1/4 case, the 
regime where 02 = vI— A?/4 gives a solution r{x) which 
is strictly positive and is the ground state of L-. Thus 
stability follows for all values of k. The last parameter 
regime, for which 02 = —1/4, is indeterminate since both 
A_ and A+ are negative and our linear stability analysis 
is inconclusive. 

These analytic predictions are confirmed in Figs, 
p^ . In Fig. ^the evolution of a dn(a;, k) solution is shown 
for 02 = 1/4 and k = 0.999. As predicted analytically, 
this parameter regime is stable for all k values. This sim- 
ulation once again illustrates the importance of an offset 
for stabilizing the condensate pOl. In contrast to this 
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stable evolution, the case 02 = —1/4 is unstable as illus- 




critical threshold the condensate destabihzes as shown for 
B = 1/2. 



FIG. 16. Unstable evolution of a Type B dn{x,k) solu- 
tion given by Eq. (^) for k = 0.5 (top panel) and k = 0.999 
(bottom panel) given a2 = — 1/4. In this case, there is no 
offset to stabihze the condensate. 



conclusively show the evolution to be unstable for all k 
values. For this case, the offset of the solution is insuffi- 
cient to stabilize the condensate. We note that for small 
values of k, the onset of instability occurs after a very 
long time. Higher values of k result in instabilities on a 
much faster time scale. Finally, we consider the param- 
eter regime for which 02 = Vl — fc^/4. In this case, the 
analytic predictions once again suggest stability for all k 
values. We do not illustrate this case since it is quali- 
tatively very similar to Fig. However, in contrast to 
the 02 — 1/4 case, for values of k close to one, there is 
a negligible amount of offset, distinguishing this stable 
case from previous ones. For these values of fc, the so- 
lution has a small amplitude compared to the potential 
so that the behavior is essentially linear and stability is 
achieved because the condensate is trapped in the wells 
of the potential, as in ordinary quantum mechanics. 
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C. Nontrivial Phase 



As stated at the beginning of Section [II, determining 
the hnear stabihty for nontrivial phase solutions is not 
amenable to analysis. This leads us to consider the sta- 
bility of nontrivial phase solutions using numerical com- 
putations. 

To begin, consider the trigonometric limit of the non- 
trivial phase solutions of Type A. These solutions are 
given by Eqs. (|l|) and (|). Figure ^ depicts the evolu- 
tion of a pair of initial conditions with Vq = 1.0 and for 
which B = 1 (top panel) and B — 1/2 (bottom panel). 
Since B determines the offset of the condensate, these 
numerical results show directly the importance of this 
offset for stability. In contrast, if the offset is too small, 
it is unable to stabilize the condensate. 

For Type B solutions, qualitatively nothing changes 
from the dynamics illustrated for the trivial phase case. 
In particular, numerical simulations can be performed us- 
ing exact solutions which are constructed subject to the 
phase quantization condition given by either Eq. (|l^ ) or 
(p3|). A numerical shooting method is used to find ap- 
propriate values of a2 for which a phase-quantized, pe- 
riodic solution exists. Once this is achieved, numerical 
simulations can easily be performed. Note that any in- 
teger value p is allowed as input for the phase quantiza- 
tion conditions, provided solutions exist for the param- 
eter values. It turns out this imposes a lower bound on 
value of p. In the simulations, the actual value of p does 
not affect the stability of the solution. Increasing the 
phase-quantization integer p leads to a solution with a 
steeper phase profile, suggesting a more unstable situa- 
tion. However, this phase effect is balanced by an in- 
creased offset a2 of the amplitude. Qualitatively, the dy- 
namics are as depicted in Figs. |K-16. Thus the nontriv- 
ial phase solutions of Type B are stable for a2 > 1/4 and 
for < 02 < Vl ~ k'^/4:, whereas the nontrivial phase 
solution is unstable for 02 < —1/4. 



IV. SUMMARY AND CONCLUSIONS 

We considered the repulsive nonlinear Schrodinger 
equation with an elliptic function potential as a model for 
a trapped, quasi-one-dimensional Bose-Einstein conden- 
sate. Two new families of periodic solutions of this equa- 
tion were found and their stability was investigated both 
analytically and numerically. Using analytical results for 
trivial phase solutions, we showed that solutions with 
sufficient offset are linearly stable. Moreover, all such 
stable solutions are deformations of the ground state of 
the linear Schrodinger equation. This is confirmed with 
extensive numerical simulations on the governing non- 
linear equation. Likewise, nontrivial phase solutions are 
stable if their density is sufficiently offset. Since we are 



modeling a Bose-Einstein condensate trapped in a stand- 
ing light wave, our results imply that a large number of 
condensed atoms is sufficient to form a stable, periodic 
condensate. Physically, this implies stability of states 
near the Thomas-Fermi limit. 

To quantify this phenomena, we consider the fc = 
limit and note that from Eqs. (|l]) and (^sj), the number of 
particles per well n is given by n = (/^ \ip{x, t)\'^dx)/TT — 
Vo/2 + B. In the context of the BEG, and for a fixed 
atomic coupling strength, this means a large number of 
condensed atoms per well n is sufficient to provide an 
offset on the order of the potential strength. This en- 
sures stabilization of the condensate. Alternatively, a 
condensate with a large enough number of atoms can be 
interpreted as a developed condensate for which the non- 
linearity acts as a stabilizing mechanism. 

Acknowledgments: We benefited greatly from dis- 
cussions with Ricardo Carretero-Gonzalez and William 
Reinhardt. The work of J. Bronski, L. D. Carr, B. Decon- 
inck, and J. N. Kutz was supported by National Science 
Foundation Grants DMS-9972869, GHE97-32919, DMS- 
0071568, and DMS-9802920 respectively. K. Promislow 
acknowledges support from NSERC-611255 



* to whom correspondence should be addressed 
[1] W. Ketterle, D. S. Sturfee, and D. M. Stamper-Kurn, in 
Proceedings of the International School of Physics "En- 
rico Fermi" (ICS Press, Amsterdam; Washington, D.C., 
1999), pp. 67-176. 
[2] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, 

Rev. Mod. Phys. 71, 463 (1999). 
[3] K. Huang, Statistical Mechanics, (John Wiley, New York, 
1963). 

[4] B. P. Anderson and M. A. Kasevich, Science 282, 1686 

(1998) . 

[5] E. W. Hagley et ai, Science 283, 1706 (1999). 
[6] M.L. Chiofalo and M.P. Tosi, Phys. Lett. A 268 406 
(2000). 

[7] Y. B. Ovchinnikov et ai, Phys. Rev. Letts. 83, 284 

(1999) . 

[8] D. Jaksch et ai, Phys. Rev. Letts. 81, 3108 (1998); 

[9] G. K. Brennen, C. M. Caves, P. S. Jessen, and I. H. 

Deutsch, Phys. Rev. Letts. 82, 1060 (1999). 
[10] D.-I. Choi and Q. Niu, Phys. Rev. Letts. 82, 2022 (1999). 
[11] G. Baym, Lectures in Quantum Mechanics, (Addison- 

Wesley, Redwood City, CA), Ch. 20. 
[12] L. P. Pitaevskii, Sov. Phys. JETP 13, 451 (1961). 
[13] E. P. Gross, Nuovo Cimento 20, 454 (1961). 
[14] D. S. Petrov, M. Holzmann, and G. V. Shylapnikov, 

Phys. Rev. Letts. 84, 2551 (2000). 
[15] D.S. P etrov, G.V. Shiyap nikov, and J.T.M. Wahaven, 

e- print |cond-mat/0006339| (2000). 
[16] L.D. Carr, C.W. Clark, and W.P Reinhardt, Phys. Rev. 



12 



A 62 0436XX-1, e-print cond-mat/9911177 (2000) 



[17] L.D. Carr, M. A. L eung, and W.P. Re inhardt. J. Phys 



[18] 
[19] 
[20] 

[21] 
[22] 

[23] 
[24] 
[25] 
[26] 
[27] 



B 33 p.xxx, e-print |cond-mat/0004287| (2000). 

M. Key et ai, Phys. Rev. Letts. 84, 1371 (2000). 

N. H. Dekker et al, Phys. Rev. Letts. 84, 1124 (2000). 

J. C. Bronski, L. Carr, B. Deconinck, and J. N. Kutz, 

PRL 

M. Kunze and et. al, Physica D 128, 273 (1999). 

Y. S. Kivshar, T. J. Ale xander, and S. K. Turitsyn, e- 



print |cond-niat/9907475| (1999). 
K. Bongs et al, 2ond-mat/0007381 (unpublished). 
M. R. Andrews et al, Science 273, 84 (1996). 
J. D. Close and W. Zhang, J. Opt. B 1, 420 (1999). 
M. R. Matthews et al, Phys. Rev. Letts. 83, 2498 (1999). 
Handbook of Mathematical Functions, edited by M. 
Abramowitz and I. A. Stegun (National Bureau of Stan- 
dards, Washington, D. C, 1964). 



[28] 
[29] 
[30] 
[31] 
[32] 

[33] 
[34] 



S. Theodorakis and E. Leontidis, J. Phys. A 30, 4835 

(1997) . 

F. Barra, P. Gaspard, and S. Rica, Phys. Rev. E 61, 5852 
(2000). 

K. Berg-S0renson and K. M0lmer, Phys. Rev. A 58, 1480 

(1998) . 



M. J. Steel and W. Zhang, e-print |cond-mat/9810284 
(1998). 

Algebro-geometric approach to nonlinear integrable equa- 
tions, E. D. Belokolos, A. I. Bobenko, V. Z. Enol'skii, 
A. R. Its and V. B. Matveev (Springer- Verlag, Berlin, 
1994). 

R. Courant and D. Hilbert, Methods of Mathematical 
Phystcs, (Wiley, New York, 1989). 

M. L Weinstein and J. B. Keller, SIAM J. Appl. Math. 
45, 200 (1985) 



13 



